I = imread('ulnar loop.bmp');

I = I(:,:,1);
I = double(I);
I = normalise(I);
%imshow(uint8(I)); hold
[m,n] = size(I);
theta = zeros(m,n);
r = zeros(m,n);
W = 7;
%[Gx,Gy] = gradient(I);
[Gx,Gy] = grad(I);
for i=(W+1)/2:m-(W+1)/2
    for j=(W+1)/2:n-(W+1)/2
        Gxy = 0; Gxx = 0; Gyy = 0;
        for h=-(W-1)/2:1:(W-1)/2
            for k=-(W-1)/2:1:(W-1)/2
                Gxy = Gxy + Gx(i+h,j+k)*Gy(i+h,j+k);
                Gxx = Gxx + Gx(i+h,j+k)*Gx(i+h,j+k);
                Gyy = Gyy + Gy(i+h,j+k)*Gy(i+h,j+k);
            end
        end
        theta(i,j) = -0.5*atan2(2*Gxy,Gxx-Gyy);
        %theta(i,j) = -pi/4+0.5*atan2(Gxx-Gyy,2*Gxy);
        r(i,j) = sqrt((Gxx-Gyy)^2 + 4*Gxy*Gxy)/(Gxx+Gyy);
        U(i,j) = r(i,j)*cos(1*theta(i,j));
        V(i,j) = r(i,j)*sin(1*theta(i,j));
        
    end
end
i = (W+1)/2:W:m-(W+1)/2;
j = (W+1)/2:W:n-(W+1)/2;
U = U(i,j);
V = V(i,j);
[i,j] = meshgrid(i,j);
 quiver(i,j,U,V)
 axis ij